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Abstract 

Traditionally, it is understood that fluctuations in the equilibrium distribution are 
not evident in thermodynamic systems of large N (the number of particles in the 
system) [la]. In this paper we examine the validity of this perception by investigating 
whether such fluctuations can in reality depend on temperature. 

Firstly, we describe fluctuations in the occupation numbers of the energy levels 
for an isolated system, using previously unknown identities that we have derived for 
the purpose, which allow us to calculate the moments of the occupation numbers. 
Then we compute analytically the probability distribution of these fluctuations. We 
show that, for every system of fixed and finite N, fluctuations about the equilibrium 
distribution do, in fact, depend on the temperature. Indeed, at higher temperatures 
the fluctuations can be so large that the system does not fully converge on the 
Maxwell-Boltzmann distribution but actually fluctuates around it. We term this 
state, where not one macrostate but a region of macrostates closely fit the underlying 
distribution, a 11 fluctuating equilibrium". Finally, we speculate on how this finding 
is applicable to networks, financial markets, and other thermodynamic-like systems. 
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1 Introduction 



The relaxation of a classical statistical system, that is isolated from its sur- 
roundings, towards the equilibrium Maxwell-Boltzmann distribution in the 
thermodynamic limit is well understood. See Huang [la] for the derivation of 
that distribution by means of Lagrange's method of undetermined multipliers 
and Friedman et al. [2] and references therein for more pedagogical derivations 
aimed at an audience that is not versed in multivariate differential calculus. 
As fluctuations in this distribution effectively all but disappear in the ther- 
modynamic limit, this distribution has proven sufficient for describing many 
thermodynamic systems for which this limit applies. 

In this paper, we consider an isolated system with a large but finite number 
of particles, N. Much work has been done on finite systems, away from the 
thermodynamic limit. For example, Hill [3 a bj has presented a wider scope for 
thermodynamics, that of describing small systems, while Lebowitz et al. [I] 
have analysed finite size effects using Taylor series expansions. More recent 
work has included topics such as the nonequivalence of thermodynamic en- 
sembles [5|6] . and phase transitions [6|Tf8f9] . 

However, in our work, we set out to derive the distribution of particles among 
energy levels/states of an isolated system with a large but finite number of 
particles, N, and the finite-A corrections to the distribution using a sim- 
ple model. Kelly [TU] has recently investigated this problem numerically, and 
macReamoinn [TT] has computed analytically the variance of the distribution 
(the amount of fluctuations). We now find closed form expressions for all the 
moments of the distribution, along with the multi-variate distribution of the 
occupation numbers. We show that at higher temperatures the amount of 
fluctuations exceeds the mean occupation number and thus no stable equilib- 
rium distribution is attained. We describe this phenomenon as a "fluctuating 
equilibrium" and we hypothesise that it could be useful to describe certain 
thermodynamic or thermodynamic-like systems, for example, in the global 
financial markets. 



2 Theoretical Analysis 

Our approach follows closely the permutational argument of Boltzmann [12] 
and is similar to methods illustrated (by examples) in textbooks [T3lll4|15j . 
However, we want more than just the most probable distribution, and we do 
not take the limit of continuity in our calculations. 

We consider an ideal gas of N particles in a volume V (y « 1) where the 
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mutual energy of interaction between the molecules is negligibly small. We 
assume that the collisions with walls of container and between molecules are 
perfectly elastic, and that the walls are rigid. For an ideal gas, the energy of 
the system, E, is given by: 

E = \Nk B T (1) 

where ks is Boltzmann's constant, and T the absolute temperature. 0. We 
assume the system is isolated from its surroundings. Following Boltzmann 
himself [12j, we assume the possible energies of the particles are given by 
equidistant, non-degenerate energy levels (with energy level spacing e), with 
the lowest energy level set to zero for simplicity. Then the energy of each 
particle can be written as a multiple of the energy level spacing, je, where 
j 6 {0, . . . M} 0. Thus we have: 

E = Me (2) 

We now use a permutational argument to describe, in a quantitative fashion, 
both the most likely distribution of the particles among the energy states and 
the fluctuations around this distribution. 

2.1 The Micro- & Macrostates 

We assume that the particles are distinguishable, as Boltzmann did. We 
describe a microstate as a certain distribution of particles among the energy 
levels where we distinguish between the different particles (e.g. by assigning 
labels to the particles) - Boltzmann called them "complexions". We define a 
macrostate as a coarse-graining of this, where we are only interested in the 
number of particles occupying each energy level. 

We denote by rij, the occupation number of the j th energy level, and specify 
the distribution of particles among the energy levels (the macrostates) by 
ft := (n 3 -)^i . The total number of particles and the total energy are conserved 



Of course, in general, temperature is normally defined as the partial derivative of 
the internal energy with respect to entropy with N and V held constant. 

2 Unlike Boltzmann, we only consider the case where the maximum number of 
energy units a particle can have is equal to the total number of units in the system, 
M. Clearly, for large N, the likelihood of a particle having energies of high j (j ~ M) 
becomes very small. 

3 Distinguishability would be a good approximation for thermodynamic systems 
away from the quantum regime of low temperature or high density. The resulting 
(Maxwell Boltzmann) statistics gives a non-extensive entropy but this does not 
concern us for the purposes of this paper. Of course, there has been plenty of 
discussion on entropy being extensive |16jl7|18"] . Also, we believe that this model 
may be suitable for economic systems where distinguishability is valid. 
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i.e. 

M M 

N = E n i and M = Y.i n i ( 3 ) 

3=1 3=0 

In this paper we are interested in the probability of macrostates. By assum- 
ing each microstate is equally likely, we can assign (unnormalized) probabil- 
ities to the macrostates by the number of microstates corresponding to one 
macrostate. This is clearly given by the following multinomial factor: 

TV! 

— ( 4 ) 

n v 

3=0 

since it is the number of ways we can place rij particles on the j th energy level 
for j = 0, . . . , M. 

The total number of microstates of the system (n(micro)) reads: 



n micro 



E » (5) 



II n 3 \ M '£^ 

N 

(6) 



3=0 

N 



y M+l \ N 



(,z . z -m I ) ( 7 ) 



2mz \ 1 — z 
1 d M fl-z M+1 \ N 



Ml dz M \ 1 - z 



2 = 



^^.cg^.iM^ (9) 

In (JSJ) we inserted the integral representation of the delta function, and in (J7J) 
and (jHl) we substituted for e~ l<i> and we computed the resulting integral using 
the Cauchy theorem. The final equality can be proven by the generalized 
Leibnitz rule and by induction in N. 



4 We have also derived a closed form expression for the total number of macrostates. 
However since we focus in this work on fluctuations around the equilibrium we will 
present that result elsewhere. 
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2.2 Fluctuations Around Equilibrium 



The most likely macrostate is obtained by maximizing the multiplicity of 
macrostates (jl]) subject to conditions ([3]). In the thermodynamic limit (N — > 
oo) this yields the Maxwell-Boltzmann exponential distribution of energies 
with the decay constant proportional to the inverse of the absolute tempera- 
ture T, as shown in [Taj, for example. For an ideal gas the absolute temperature 
is proportional to the mean kinetic energy of the particles which, in our model, 
is proportional to the ratio jr. This can be seen clearly from equations ([Tj) 
and (J2J). In fact, we have: 



f = ffT (10) 



In what follows, we use "temperature" and "T" for this ratio . Of course, 
this ratio is the "specific energy", simply the average energy per particle in 
multiples of the energy levels spacing e, which is a well defined quantity in 
itself, no matter how we define temperature. 



While the most likely distribution is well known, little is known about the 
fluctuations around this distribution. This raises interesting questions. The 
traditional understanding of statistical mechanics suggests to us that an iso- 
lated thermodynamic system will always move to the Maxwell-Distribution as 
a result of the Law of Large Numbers (LLN). We wish to investigate the effect 
of the temperature (the ratio ^) on this convergence. If the energy content of 
the system is "very large" (see later for a comment on what this might mean), 
such that the temperature is high, maybe the best the system can achieve is 
only an equilibrium region of states; "a fluctuating equilibrium". Indeed, at 
very high temperatures, the fluctuations could be so dominant that disorder 
and chaos would prevail at the macroscopic level. We believe such behaviour 
may be apparent in thermodynamic-like systems where N is much smaller 
than Avagadro's number. 



In order to investigate this possibility we must first compute all the moments 
of the occupation numbers of the energy levels. We have: 
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(11) 



p=0 



dlog(x) 



cf' 



Jge^ J2 l^fi( c "^(M^-l) + l))(12) 



n n ' p =1 
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rflog(x) m 7 2?r 



p=0 



e ^ 1 (e J ^. 



A? 



\ 1 - e - '? 
dm ! d u (i_ z m+i 
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d\og(x) m M\ dz M \ 1 - z 
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x=l 
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AT-1 



rflog(x) 



X: i 9J <m(x - iyc^c^:r-T + <W(* - 1) 



q=0 



x=l 



(N-l)Am 



(m) 
A/ 



x=l 

(13) 
(14) 
(15) 
(16) 



In ( TT2|) we introduced an integral representation of the delta function and we 
used the identity: 



cf 



n" 



d\og(x) 



(17) 



x=l 



for m a non-negative integer. In (flBl we computed the sum using the multino- 
mial expansion formula and in (|T4"|) we substituted for z = e~ l<?i and computed 
the resulting complex integral using the Cauchy theorem. In (1151) we computed 
the derivative using identity (1A.8|) and in (FIB"!) we used the identity (lA.lOj) . 
The coefficients 4 m) are defined in ([Q2l> - (j04j> . 



From (Ti~6"l) we see that for iV — > oo the moments of the occupation densities 
of the energy levels, Xj := rij/N, read: 



•i' ; 



M-mj+N-l-m 



lmj<Af 



°AT-1 



M+N-l 



°iV-l 



N' 



and are all finite, once the temperature T = M/N, is finite. 



Now imagine that N — > oo but T is kept finite. Then from (1181) we have: 
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M + N 



inj 



[1 + T)"- 
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UtTT] 



-imlog(l+T- 1 ) 



(19) 



and we conclude that the moments of the distribution depend on both the 
energy of the level j and on the temperature of the system. In the case m = 1 
we retrieve results from [2] (equation (8), page 119) and from [10] (equation 
(6), page 10). In the case m = 2 the result fits in with [TT] (equation (33), 
page 9). 



Thus all the moments, and in turn also the distributions of the occupation 
densities, depend on the temperature T only, and not on either iV or M alone. 
In order to check that the formula is correct we use MATHEMATICA to show 
that the total number of particles is conserved: 



M 

E 

3=0 



\ x i) 



E ' V ~ 2 
3=0 



+N-2 



U N-1 



N-l 



(20) 



For the limit N — > oo with fixed T (which means that M — > oo as N = TM) 
we get: 



M i M 

E fa) = TfTrY S e-^ l+T ^ 

3=0 \ L -r 1 ) 3=0 

1 1 ~ (ijf 
(1 + T) 1-^ 



as expected. 

To quantify the fluctuations we look at the variances of the occupation densi- 
ties, which from (fl6l) read: 
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U N-2 / 1 L/ N-2 



L j<M ^jM+N-l ^jy qM+N-1 



(N - 1) C ™- 2j+N - 



+ l 2j<M ^7 „M+N-1 ( 22 ) 



1 1 T j I T j 



Equation (122]) follows from the definition of the variance and equation (1231) 
follows from ( fT8]) and ( fl9l) . For given iV and j, the variance is bounded from 
above as a function of T and its maximal value, which occurs at the temper- 
ature T = j, is given by: 



[af^f = i (of ax ) (1 - (af ax )) (24) 



for ^a;™ 3,5 ^ = + Note that, for j > 0, both the mean occupation 

density and the variance tend to zero at high and low temperatures, yet at 
different rates. Therefore we consider the ratio -Ar. We get: 



(25) 

(xj) (xj) I i T \ if t^oo 

\ V 

Thus, in the limit of high and low temperature^, the fluctuations prevail, the 
system diverges and the equilibrium distribution does not exist, contrary to 
the traditional view of classical statistical mechanics [lb] . We plot the stan- 
dard deviation of the occupation densities in units of their mean occupation 
densities as a function of the temperature for different values of N and j 
in Figure IF. 1L In addition, we plot the occupation density and its standard 
deviation for j = 1 separately in Figure IF.21 

We conclude this section by computing the distributions of the occupation 
densities, meaning we compute the likelihood that the number of particles Nj 
at the j th level equals n,-. We get: 
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As T — > both (xj) and <jj approach zero for j > 0, but (xj) does so at a faster 
rate. 
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N-l 



22=0 



l qj<M^hA- 1 ) ^q^N-l-q +°Nj,M^nA~ l ) 



(26) 

(27) 

(28) 
(29) 



In the first equality in (|2"TI) we repeated the manipulations from mean- 
ing we have inserted the integral representation of delta functions, performed 
the sum using the multinomial expansion formula and finally evaluated the 
complex integrals over a pair of unit circles using the residue formula. In fl28l) 
we have computed the derivative by Z\ using equation (IA.8[) and in fl2"9"|) we 
computed the derivative by z 2 using derivatives of elementary functions. 

For j < T the large N limit of the distributions of the occupation densities 
reads: 



N 

Jim P (iV, = nj) = £ C^C-ir^Cf fo) 4 



g=nj 

Cg foP (1 - (30) 
Normal(A^ (acj) , iV (a;,) (1 - (xj)))(n) 

/ (n iV X 



Equation (|30|) follows from the fact that the large A" limit of the last binomial 
factor on the right-hand side in (|29|) reads (xj) q and equation ([Mil follows 
from the normal approximation to the binomial distribution. We plot both 
the exact result and the large A" limit in Figures IF. 31 and IF.4I for several values 
of A" and M. 

Thus the distributions of the occupation densities, Xj = rij/N, conform to 



Gaussians with means A" (xj) and standard deviations l/y/N y (xj) (1 — (xj)) 
in accordance with the results from [la], for example. However, for every finite 
A" if the temperature T is high enough the occupation number in units of the 
mean occupation number has a distribution whose "flatness" or "half width' is 
arbitrarily high. In other words, for every finite A" the amount of fluctuations 
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diverges when the temperature gets very large. We depict that statement in 
Figure IF. 51 



2.2.1 A Note on Units 

In this paper, we call the ratio ^ the "temperature" and discuss how high 
values of this quantity can lead to more fluctuations in the distribution. Equa- 
tion (fTUj) relates this ratio to the absolute temperature and the "energy level 
spacing" , e, for the ideal gas system. But when is this ratio "high" ? That is, at 
what absolute temperature does this occur? Suppose we take the case where 
M = lOOOiV, and we consider e as some sort of precision limit on measure- 
ment e.g. ~ O.leV ~ 1.6 x 10~ 20 J. Note that this is the amount of energy 
in an infrared light quantum and, as such, is measurable by spectroscopical 
methods. From fflQj) this gives an absolute temperature of about 770,000 K 
which is certainly hotter than anything naturally occurring on Earth, though 
would be easily exceeded in most stars. Thus our analysis would not be use- 
ful for physical systems at room temperature although may be of interest for 
astrophysics. 

However, suppose now we consider a financial system, and let the "energy 
quantum" represent now a unit of money, e.g. one dollar (with the Boltzmann 
constant being suitably redefined or discarded). A typical financial system, 
where the number of market participants is much less than the amount of 
money involved, could give values for the ratio of 10 4 , or much higher. The 
considerations in our paper would therefore take on greater importance for 
these systems. 

2.2.2 Multi-point (Macrostate) Probability Distribution Function 

We now generalize the calculations in (r2"6"j) - f[2T?j) to obtain the multivariate 
probability function. We take p = 1, . . . , M + 1 and an ascending integer 
sequence < j\ < ■ ■ ■ < j p < M and we define the p-variate probability 



function P I f] (Nj s = nj s ) I as the probability that we find hj s particles at 



energy level j s for s — 1, . . . , p. The function is given in fl32l) . We have: 





N ~ 1 / P (q) t ^ 

q=0 \l = l 

( V (jv) 

\i=i 




m 




) 



(JV) 



, Nr ,M-J*+N-l-q 



) 



(32) 
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subject to: 

q N 

m r 9 ^ = ^2 ^, S; for r = 1, . . . , p and m^' = <5 rjSi for r = 1, . . . , p. 
i=i i=i 



Here, we have defined: 



J' 



(33) 



Note that for every q = 0, . . . , N — 1 the expression on the right-hand side 
is summed over integer grid-points (sj)f=i °f a (/-dimensional hypercube of 
side length p i.e. (sj)?=i £ {1, 2, . . .p}. As such the expressions in parentheses 



1 )i=i 

These numbers count how many coordinates of the grid-point are equal to 

I = 1, . . . , q. Therefore the p-variate probability function is invariant under 

p , q \ 

permutations of its arguments. In addition, the identity m i — Q holds 

i=i 

true. The proof is in Appendix C 



contain p q terms. The expression on the right-hand side depends on (m 



In the case p = 1 we have S\ = ■ ■ ■ = = 1 and mj = iV and j Sl = ■ ■ ■ = 
js q — ji — 3i the expressions in parentheses reduce to one term only, and we 
clearly retrieve the probability function in (1291) . 

Expression (132]) is cumbersome and does not provide much insight in to the 
problem. Thus, in order to better understand the expression, we compute its 
large- N limit. We prove that the p-variate probability function is normalized 
to unity in Appendix E In addition we take j s = s — 1 for s — 1, . . . ,p., i.e. 
we are interested in the p lowest energy levels. The result is a multinomial 
distribution with likelihoods of individual trials pi given in (|35l) . We have: 



p n w 



Js — n js, 



\S = 1 



p+1 

Up, 



(AT- fi)!(n» A '=i 

1=1 



(34) 



where 



T \l 1 



Pi 



i) t in = i,...,p 



5$j)f if I = p + 1 



(35) 



The proof is in Appendix D 



2.2.2.1 Specific cases 

In the case p = 1 the expression (ID. 61) reduces to the binomial distribution 
with mean N (xq) and variance N {xq) (1 — (xq)) in accordance with equation 
(!30l) . In the case p = M + 1, from (1D.6[) . (1351) and the conservation laws (J3]) 
and by using the Stirling approximation to the binomial coefficient, we get: 



11 



M 



Here we have used um+i = N — ^ n» = 0. The last line in ( 136|) is, except for 



a multiplicative prefactor, equal to the number of microstates corresponding 
to one macrostate, divided by the total number of macrostates (compare (HI) 
and ([9])), which is what we would expect. The fact that the prefactor is not 
unity is due to the approximations made (large N approximations, Stirling 
approximation etc.). This could be improved upon using higher orders but is 
not needed here. 



2.2.2.2 The mean, the covariance and the correlation matrices 

Now take p — M + 1. In order to quantify the entire amount of fluctuations 
we give the mean and the covariance matrix of the (M + l)-point distribution 
function. The mean reads: 



in accordance with the mean of the univariate distribution, given in (Tl9l) . The 
covariance matrix reads: 



for li, l 2 = 0, . . . , M. in accordance with (123]) . The results (1371) and (!38|) follow 
from the fact that the multivariate distribution is a multinomial distribution. 

From (1371) and (|38|) we compute the Pearson correlation coefficient be- 
tween occupation numbers at levels i ^ j. 



(m) = Np t+1 



N (xi) for I = 0, . . . , 



M 



(37) 




(38) 
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Pm,rij 



{{rii - (n,))(n j - (n,))) 



^(I-^XI -<*.)) 

_y(i+i)/2 for T -»• and i,j > 
— T _1 for T — > oo 



(39) 



We can see that occupation numbers of any two different levels are anticor- 
related. The occupation numbers become uncorrelated once temperatures be- 
come high or low and a maximal correlation is attained at some intermediate 
temperature (see Figure lR6l) . This is what one would intuitively expect, since 
at both high and low temperatures the mean occupations tend to zero and as 
such can be treated as independent from each other. 



2.2.2.3 "The total amount of fluctuations" in units of the mean 
occupation. 

We conclude this section by investigating how strong is the effect of fluctu- 
ations on the whole vector of occupation numbers rather than only on one 
particular occupation number, which was the problem that we explored in 
(1251) . Here we need to use some measures associated with the covariance ma- 
trix and the vector of mean occupation numbers. We choose to compare the 
square root of the trace of the covariance matrix to the L 1 norm of the vec- 
tor of mean occupation numbers. Note that this is a natural generalisation of 

(\\y-(y)\\ 2 ) 

the variance for multivariate distributions, i.e. it is of the form: — ^4y| — 

where |.| and |.| are the L 2 and L l norms respectively and Y is the random 
vector. We define g := (cij)f^ =0 and n := ((rij))^ , and make the substitution 
j: = jyrpj-. Using (l3"Tj) and fl38l) . we calculate the following: 



t(g) j ~ x / 2 - ? A/ _ p A/+l + J.2M+1 _ J.2A/+2 



\n\ V( 1 + f) ViV(l-r M+1 



N 



-r± if T -> oo 

V±-e~ N (41) 
VT if T -> 



We examine the approach to the high temperature limit more closely by graph- 
ing (in figure IF. T[) the function in equation (j4"0l) against T, for different values 
of N. We can see that, for fixed N, the total amount of fluctuations increase 
with T towards the constant value given in equation (j4"Ii) . as opposed to the 
univariate case in equation (1251) where the "amount of fluctuations" diverged 
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as T — > oo. Therefore, the conclusion from (|25|) still holds and the system 
does not converge on the Maxwell-Boltzmann distribution but instead fluc- 
tuates around it except that the total amount of fluctuations asymptotically 
approaches a finite limit. 



2. 2. 3 Corrections to Fluctuations for Finite N 

We can see, from fl3Tl) . that in the thermodynamic limit, N — ► oo, the amount 
of fluctuations present in the system depends on the temperature T and the 
total number of particles N. For finite systems, however, there are corrections 
to the amount of fluctuations which we specify now. We note that the last 
term in (1291) is absorbed as the q = N term into the sum and the sum over q is 
performed by expressing the ratio of binomial factors as a Laplace transform 
by q viz: 

Ui<M N -'m%-i = =: / d9u j (9)9« (42) 
u v-i J 

where the weight 0Jj{9) reads: 

M+N-l-q(j+l) 

»m = J y^w (43) 



Exchanging the integration in ( )42l) with the sum over q in (1291) we easily arrive 
at the following result: 

p (Nj = hj) = f 1 dBuJ d {Q)e n i (1 - e) N ~ n ' (44) 
3 Jo 



The integral in (T4*3|) can be done analytically by using complex calculus, how- 
ever, since we see in Figures IF.3I and IF. 41 that the corrections are not higher 
than a couple of percent even for iV of the order of a couple of tens we do not 
see any motivation in deriving these results analytically. Thus, we conclude, 
the distribution of the number of particles on the j level is a continuous linear 
superposition of Binomial distributions with mean NO and variance N6(l —6) 
and weights u)j(9) for 9 £ [0, 1]. In the thermodynamic limit the weight tends 
towards a Dirac delta function that picks out 9 = (xj) and we retrieve the 
result from equation ( l30i) . 



3 Conclusions 



It has long been understood that the state to which an isolated thermody- 
namic system is attracted to, as a result of the Law of Large Numbers (i.e. 
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N — > oo), becomes an exponential distribution in the thermodynamic limit. In 
this document however, we have shown that for every finite N (large enough 
for some of the formulae simplifications used in this document to be valid), 
there always exists some degree of fluctuations that depend on the tempera- 
ture of the system. The occupations of the energy levels fluctuate, and their 
distributions can be well approximated by normal curves with means N (xj) 
and variances N (xj) (1 — (xj)) (here (xj) is the mean occupation density of 
the energy level and is given in (fl8l) ). In fact for higher temperatures the 
fluctuations can be considerable relative to the mean occupation, as seen in 
Figure (IF.ip and Figure (IF. 21) . Thus, under these conditions, one cannot speak 
about a static distribution of particles among the energy levels since in any 
finite fixed N system, there is always some degree of fluctuation about the 
underlying mean distribution. 



In this work we used "temperature" in the sense of a specific energy (energy 
per particle). If we look at other systems, this temperature could be a proxy for 
any "mean quantity" (i.e. the amount per component, of a conserved quantity 
that is distributed among the components) in any finite N thermodynamic- 
like system. This understanding that there are increasing degrees of disorder 
associated with an increasing mean quantity can be applied to any average 
quantity. 



Furthermore, several people [19 2 0l[2Tl[22] have recently applied Boltzmann 
statistics to income distributions and the financial markets. We believe our 
work can help explain why financial markets are an economic equilibrium that 
is not static, but fluctuating. We believe this work has the potential to lead 
to a method for measuring the temperature (i.e. the amount of fluctuations) 
and the rate of change of temperature (i.e. trending) in financial and property 
markets. 



But possibly the most interesting thing of all is to think about what happens 
in a thermodynamic-like system, when the components of the system are not 
simple things like particles, but entities with some intelligence having the 
ability to adapt to their environment. Such adaptive entities would be capable 
of choosing to convert potential energy stores into kinetic energy (and vice 
versa) and as a result the overall system would effectively have the ability to 
"endogenousltf alter its own temperature. 



It has long been understood that simple systems spontaneously move to a 
thermal equilibrium. In future work we hope to show that (as a result of 
competition and cooperation) " complex adaptive systems" gravitate to another 
type of equilibrium, an equilibrium somewhere between order and disorder, an 
equilibrium which maximizes survival by maximizing flexibility, an equilibrium 
which we might call a " structural 1 or " organizational equilibrium" . 
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Appendix A 



Here we list some important identities that are used in this paper. We believe 
that these identities are useful per se in analysis, partial differential equations, 
and probability & statistics. 



The Sum of Powers Identity: 



t-1 n+l 

£* n = £ 

1=0 k=0 



£'"-E<"-' +1 7^W* (Ai: 



where 



k 



1 1 



c fe :=£(-l) P £ ~ = ( 1 '-2'T2'-R'-) (A ' 2) 

p =o ni+ ... +np=p+k i\ n q \ v 

ni,...,n p >2 q=l 



for k = 0, . . . , n + 1. The identity ([Oil can be proven by using a trick 
dfog^nt 1 t = l n for n, I G N, changing the order of differentiation and summa- 
tion, re-summing the resulting geometric series, and differentiating the result 
using the chain rule of differentiation. 



The Combinatorial Identity: 



E II ^K-T-l = C t S n-s (A.3) 

l<h<-<l B <4-n j=a 

with l = p = and l s+ i = p s+ i = 5. The identity (1A.3|) is easily proven from 
combinatorial considerations. The proof is left to the reader. 

The "Sum Over a Simplex I" Identity: 
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e mr = ,"'- a - irw ' (a.4) 

a<l s <-<h<b 1=1 II (i + (-W s - iV"- J )) 



=1 



where := £ ra^ for j = 1, . . . , s. The sum (IA.4I) is done by proceeding 

9=1 

from indices with large subscripts towards those with small subscripts and at 
every time applying the identity ( 1A.1I) . In doing this we retain the highest 
order term only. However, inclusion of lower terms is possible since they are 
of the same form as the highest order term. It is only that the enumeration of 
all possible terms that emerge is cumbersome. This is left for future work. 

The "Sum Over a Simplex II" Identity: 



p-i 



0<no<---<np^i<np q=0 

P yj n P +p-l 
p 11 a l 

= U-lY-i (A.5) 

9=0 n (1 - ai a g _0 n (1 - a q a t ) 

1=0 l=q 

The identity follows from a iterative application of the geometric sum formula. 
The "Integral Over a Simplex" Identity: 



0<&<-<£i<l j ~° 

n mf\ ■ n f ^ \ (A.6) 

1=0 1=1 [j - + E rnq 1] + mf 

\ g=i-j j 

where = and £o = 1- The integral (1A.6I) is computed by integrating 
in decreasing order of the subscript, i.e. starting from £j and ending at £i, at 
each time substituting for £j = where tj G [0, 1] for j — i, . . . , 1 and by 

using the identity : 

dtt n {\ - t) m = - — — (A.7) 

v ' (n + m + 1)! K ' 

o v ' 

for n,m > —1. 

The "Power of the Sum" Identity: 
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M I N-l \ 

T, zP [(H W**af^) + <W* + (A.8) 

p=0 \ q=0 J 

where the coefficients read: 

MN) = QNgr-u+N-l-q (A>g) 



for q — 0, . . . , N — 1. The proof is in |Appendix B 
The "Differential" Identity: 



's-l 



Jin m 

' V - If = E ^ - ir s e sy I I fl i m) (A.10) 



,=i v=o 
where the coefficients satisfy recursion relations: 



l (™+ 1 ) = sa^l s < m + aKl s > 2 

= e nV-o n! (a.ii) 

s-l «=o 
n ; =m— s+1 

i=0 

E (l) m " s+1 11(^737)" ( A - 12 ) 

0<no<-<n«_ a <m-s+l Z=0 S 

8-2 

n/ s— I \m— I 
, \s-l~\> 

s-l 



= t^W EM)^""^ " <?) m ( A - 14 ) 

l S J-J- g=0 

= {{1}, {1,1}, {1,3,1}, {1,7, 6,1}, {1,15,25, 10,1}, 

{1,31,90,65,15,1}...} (A.15) 

with = 1. The equality in (1A.11I) follows from iterating the recursion re- 
lation and recognizing a pattern in the consecutive iterates. The equality in 
( 1A.12I) follows from parameterizing the sum over the simplex and the equal- 
ity in (IA.13P follows from flA.5[) . Finally, the equality in (IA.14I) results from 
simplifying expression (IA.5|) . In flA. 15[) we give numerical values for the co- 
efficients for m = 0, . . . , 5. Note that since recursion relations such as (1A.11I) 
appear in such vast fields as anomalous diffusion processes and econometrics 
& time series modelling, mathematical techniques for solving these relations, 
presented here, are valuable for researchers from those fields. 
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The "Polynomial" Identity: 



n n 2p 

(«+i)„- 1 :=n(«+j) = E« n - p ^r 

3=1 P =0 Z P- 



(A.16) 



The identity fl A. 16[) is valid in the limit n — > oo and it follows from expanding 
the product into a sum and computing the coefficients at powers of a using 
identity fDOj) . 



An Auxiliary Lemma: 



n-p2 



m=pi 



(n-p 2 ) a2+1 n ai 



J n 



J n 



(A.17) 



n 



rPl 



, a 2 + 1, ai + 1) - a 2 + 1, a x + 1) (A.18) 

n n 



Here 5 is the incomplete Beta function. The result flA.17|) follows from ex- 
panding the term in parentheses in a binomial expansion, doing the sum over 
ri\ using the leading order term in flA.lj) and then re-summing the binomial 
expansion. 



The Generalization of the "Auxiliary Lemma": 



N 



°n-p 1 ,...,p N 



niH hrtjv=n j=l 



n 9=1 



d N ~ x e 



JV-2 2^ a i 

HOf 1 (l-^-) aw_i+1 



V 



'JV-1 



^(l-^-ir n 1 



J=1 ""j-i - J - n0 j-i 



The result ( 1A.19I) follows from an iterative application of ( 1A.17I) and perform- 
ing the manipulations that were used to derive ( 1A.17I) . Note that when all 
Pj = the integral on the right hand side factorizes and equals the multivari- 



N 



ate Beta function, which is Yl 
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Appendix B 

We prove the "Power of the Sum" identity (IA.8I) : 



1 - z M 



M+l \N ( M N V 



1 - Z 

M p 

E- P X 

\p=0 pi=0 



■u) 



\ N-l 

£^ E(W*+ 1 )(W,*«+ 1 )) = (■■■(5p,2 j u 2 + 2-l j < p u+p+l)) N ~ 1 (B.l) 



AT-2 

E (<w« + i)(^- P1 ,2i« 2 + h< P - P1 2u + cr pi+i ) x 

Pi=0 

(5 P ,3i^ 3 + 3 • l 2j < P u 2 + 3 • l^CP^ + C P 2 +2 )) N ^ (B.2) 

\ JV-4 

£ (S Puj u + l)(5 p - pu3j u 3 + 3 ■ l 2j <p- Pl u 2 + 3 ■ l J < p . Pl uC p l - pi ' J+1 + CT P1+2 ) 

Pi=0 

(■ ■ • (5 pAj u* + 4 • l^ p u 3 + 6 • l 2l < p Cr 23+1 u 2 + 4 • l^CT^ 2 + Cf+ 3 ))^ 4 (B.3) 



Thus identity (1A.8|) holds for A" = 2,3,4. Assume identity is valid for any 



value of N. Then the coefficient at z p in S N+l reads: 



P N-l \ 

e + 1) e ^< P - P1 4™ ,AV + <w*i«* (b.4) 

pi=0 \<Z=0 / 

'iV-1 \ (N-lp-qj \ 

E i (9+ i),-<^^^)^ +1 + + E E 4™<^V + i^u" 

? =0 / \ ? =0|ii=0 / 

N p-qj 

2^Ugj< p O (/ _ 1 L- Ar „ (2 lg>i + l 9< Ar 2^ % ^JV-l-g J M + ±Nj<pV> + <J(N+l)j,pU 

q=0 pi=0 
N 

X^(1 r N nP-lj+N-l} _i_ 1 r<Nfyp-qj+N-q\ q . 1 JV , r 

q=0 

N 

g=0 

AT 

E^'^^' + V);/^ q- e - d - (B-5) 

g=0 
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Appendix C 



Here we derive a formula for the p-variate probability function of the occu- 
pation numbers, i.e. the likelihood to find n Ja particles on energy level j s for 
s = 1, . . . ,p, respectively. The calculations are a natural extension of calcula- 
tions (l2l)]) - ([29|) . in the p = 1 case and thus we leave the explanations to the 
reader. We take an ordered sequence < j% < ■ • • < j p and we write: 



p 
1 

s=l 



E 



M 



M 



n n ' M >x«* V»=i 



Tl 6 n ja ,f 



1 d 



q=0 

W / /' i ,//< \ / i ^M+l p 



1 d n ^ \ l-zf 



N 



0s 



v 1 d**« 



1 - Z l s=l 
N-l <? 



Q>i*(2i +s -i)) 



21=0 



n 



b 9 J-jf^M^jV-l 



M-J?+N-l-q 



(C.l) 



(C.2) 



n 

i=i 



Zl + s=0 

N 

i) + w -(^d 



2=1 



(7=0 



\Z=1 



(C.4) 



N 



subject to = J2 5 rSl and = 5 r Sl for r = 1, . . . , p. Here J„ '■= J2 js t ■ 

i=i i=i ' H i=i 



Appendix D 



Here we prove the large- N limit of the multi-point probability function P I f] (Nj 



s=l 



for the case of j s — s — 1. Recall that (jz)f =1 is a strictly ascending sequence 

_ q 

< ji < ■ ■ ■ < j p < M. Then, from (1331) . J% = (J2 si) — q and we have: 

i=i 
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p p r M-J s +N-l-q 



M+JV-1 



«l=l s q =l 

# of r's in (S;)' =1 = rn r 
for r — 1 , . . . , p 



/"fin -H 
U AT-1 



<3 

p p 1 / T \2> ! 

# of r's in («0?=i = ™r 
for r — 1 , . . . , p 



(D.l) 



1 q\ { T \g> 



(D.2) 



i=i 



In ( ID. If) we used ( |T8l) and ( |T9l) . For simplicity we drop the superscript in the 
m indices from now on. Note that the sum runs over all integer grid-points of 
a g-dimensional hypercube of side length p subject to the gridpoint having m r 

coordinates equal to r = 1, . . . ,p. In ( ID. 21) we parametrized (sz)Li := U)j=i,i=i 

r 

for r = 1, . . . , q subject to J2 Po — rn 9 for 6 = 1, . . . ,p with p l e > and 

i=i 

q p 

we summed the p parameters. Note that since the exponent si = J2 Imi 

i=i i=i 

does not depend on the p parameters the term in sum is multiplied by the 
cardinality of the set to be summed over and the cardinality in question equals 
the multinomial factor. Thus, from (|32|) and from (ID. 21) . we have: 




In ( ID. 41) we defined 



n 



nj t and in (ID. 51) we performed the sum over 



the rrii parameters using the multinomial expansion formula. In (1D.5[) we 
performed the sum over q using the binomial expansion formula. Thus, as seen 
in (1D.6[) . the result is a multinomial distribution with likelihoods of individual 



i=i 
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trials given as: 



, if i= i,..., P 



We check that the likelihoods of the individual trials sum up to unity. We 
have: 



2-^1 rl rp I I 1 

i=l J \/ = l J 1 



t + r 



l-( ) p + ( ) p = l (D. 



as expected. 



Appendix E 



Here we prove that the multi-point probability function (132]) is normalized 
to unity. Firstly, for given q we parametrize the sequence of s's as follows 



v 



i 



i s i)i=i := S 2 &i k° ) • This means that we assume that the sequence of 
\e=lpg=l * p °J l=1 

s's contains Ti occurrences of the integer (i = {1,2, . . . ,p}) at positions 

v 

k\, . . . , fe~ . Hence J„ = J2 r eje and ml = r t for I = 1, . . . ,p and ml = 

<?=i 

otherwise. Inserting this into (|32|) and summing over the occupation numbers 
we obtain: 
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( 



n 



N jp = n jp ) 



N-l 



E E 

5=0 P 

E r 8=1 
8=1 

ri>0,...,r„>0 



V s=1 e=i 



^-(E r e i e )+Ar-l- 9 



+ E 5p 

E r e =Ar 9-1 0=1 



A ' niEcS3.(-i) r '-^'^.< ( , 



0=1 
n>0,...,r p >0 



N-l 



( 



E E 

5=0 p 

E r «=9 

n>0,...,r„>0 



1 p 

0=1 



\ 

p Af~(E r e j e )+N-l-q | 

ri(i-ir e -cf^-i-T 



! 0=1 



/ 



AH X 

• na-i) 



+ E ^ 

P E ff r J "l 



n>0,...,r p >0 
vM+iV-1 



M+iV-1 



J-0<M^Ar-i + fo,M"o,Af — ^N-l 

In ( 1E.2fl we summed over the individual occupation numbers using the bi- 
nomial expansion formula and the fact the term on the right hand side fac- 
torizes. It is readily seen that out of the entire sum in (IE. 21) only the term 
at ri = • • • = r p = along with q = and N = in the first and second 
expression respectively is picked out which leads in a straightforward way to 
(IE.20 . This finishes the proof. 



(E.1) 



(E.2) 



(E.3) 
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Fig. F.l. The standard deviation of the occupation numbers of the energy lev- 
els in units of the mean occupation numbers as a function of temperature T for 
N = 10, ... , 10 5 and j = 0, . . . , 5 (with growing dash length). At low and high 
temperatures the fluctuations are of the order of a couple of units of the mean 
occupation, thus "the system diverges" . 



( j-l, N-10) ( j = l, N=100) 
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Fig. F.2. The occupation density (thin line) of the first energy level, along with its 
standard deviation (thick line) as a function of temperature. Here the number of 
particles reads N = 10 (left) and TV = 100 (right). We see that at low and at high 
temperatures the standard deviation exceeds the mean occupation density and thus 
the system fluctuates strongly. 
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Fig. F.3. The probability distributions, P(rij), of the occupations numbers of energy 
levels j for both the exact formula (thin line) and the formula derived in the ther- 
modynamic limit (thick line). Here j = 0, . . . , 3 (from right to the left) and N = 50 
and the value of T in the plots clockwise from top left is T := M/N = 1,2,3,4 
respectively. 
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Fig. F.4. Same as in Figure lK3l but now have iV = 100. 
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Fig. F.5. The probability distribution of the occupation number of level j = 1 for 
temperatures T = 10, 20, 50, 100 and for system sizes TV = 16, 64, 256, 1024. For 
every finite N the system "diverges" when the temperature goes to infinity. 
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Fig. F.6. The modulus of correlations between occupation numbers at different levels 
i = 0, . . . , 5 (with increasing dash length) and j for j = 1, . . . , 4. Occupations of 
two different levels can be roughly treated as independent from each other unless 
i,j = 0, 1 or i,j = 1, and temperatures are close to unity. 
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Fig. F.7. The total amount of fluctuations of the vector of the occupation numbers, 
graphed as a function of T, for N = 10, 30, 50, 70, 90. The amount of fluctuations 
were expressed as the trace of the covariance matrix of the occupation numbers in 
units of the L 1 norm of the vector of mean occupation numbers. 
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